Tutorial: nf-selection

Introduction

In 1858, scientists Charles Darwin and Alfred Russell Wallace described for the first time and independently, the principle of natural selection [1]. This revolutionary theory, published in the book “On the Origin of Species by Means of Natural Selection”, articulates the idea that beneficial traits — those characteristics that improve an individual’s chances to survive and reproduce — tend to become more frequent in populations over time [2,3]. A few years later, the Darwinian concept of natural selection was combined with the rediscovery of Mendelian laws to define selection as a change in the allele frequencies of a population by a specific selective pressure [4,5]. The redefinition of the concept of selection from a genetic perspective gave way to the possibility of studying and identifying the impact of selection at the most fundamental level, that is, the genome [4].

Over the past few years, the considerable increase in computational power and the development of next-generation sequencing technologies have given way to the opportunity to identify characteristic patterns of selection in the genetic material of various species and organisms [6]. In doing so, researchers have begun to understand not only the evolution of the populations but the geographical distribution of traits and diseases as well [7,8,9]. Although detection of selection signatures increases our understanding of populations, the required use of statistical and bioinformatics methods represents a barrier for scientists unfamiliar with data manipulation. To overcome this problem, several bioinformatics pipelines have been developed to simplify the steps required to move from genotype data to hypothesis regarding signatures of selection. For instance, the web-based tool “Selection Browser 1.0” has the purpose of determining the action of selection in data from the 1000 genomes project by implementing a series of informative statistics 10. Likewise, Murray Cadzow and collaborators created “Selection Tools”, a pipeline to take Variant Call Format (VCF) through several selection analyses [8]. Although these resources allow the simplified treatment of the genotype data, both are based on the same detection methods - Tajima’s D, CLR, Fay and Wu’s H, XPEHH, iHH, iHS, FST, DAF, and XPCLR -, thus leaving the global overview incomplete.

The Population Branch Statistic (PBS) is a population differentiation-based method with strong power to detect recent natural selection [11]. The Integrated Haplotype Score (iHS) is a linkage disequilibrium-based method to detect recent positive selection. Both statistics have been used together for detecting regions of the genome under selection in several investigations [See references 12-15]. Despite this, there are currently no public bioinformatics tools that allow PBS and iHS to be calculated in an automated way. Here we provide a brief overview and tutorial for nf-selection, a simple Nextflow pipeline for identifying recent signatures of selection in whole-genome genotype data with PBS and iHS. This pipeline compute PBS for any diploid species where genome-wide data is available in a VCF file.

Methods for detecting signatures of selection

The methods for the study of selection can be divided according to the temporal depth over which the selection occurred, the type of selection being investigated, and the approach of each of the methods. This last classification includes frequency-based methods, linkage disequilibrium-based methods, and population differentiation-based methods [3].

Population differentiation–based methods have been widely used to identify selection signals on a timescale in the range of 50,000 to 75,000 years from the present [8,16]. The foundation behind these statistics lies in the fact that separate populations experience different selective pressures, and, as a result, the adaptive traits of one population may be dissimilar from the traits of another. Thus, if the selection is acting at a specific locus of a population, but not at the same locus in a different population, we will observe significantly different allele frequencies. Unlike other approaches, population differentiation–based methods have the ability to detect positive, neutral, and negative selection [3].

The Wright’s Fixation Index (FST) is the most widely used population comparison-based statistic for the detection of selection signals [16]. This metric is responsible for comparing the variation of allele frequencies within and between populations. If there is a marked differentiation at one locus between two populations, FST values will be relatively high. By contrast, if two populations are homogeneous at one locus, the FST values will be relatively small [3]. Over the last few years, a series of statistics have been derived from the FST, among them, we can find the Lewontin-Krakauer Test (LKT), the Locus-specific Length Branch Metric (LSBL), the Cross-Population Composite Likelihood Ratio (XP-CLR), and the Population Branch Statistic (PBS) [3].

Linkage disequilibrium-based methods detect selection signals by searching for regions with long haplotypes. This type of statistics are particularly useful for detecting variants that have undergone partial or incomplete selective sweeps (variants that have not been fixed in the population). A large percentage of LD-based tests derive from the Extended Haplotype Homozygosity Statistic (EHH). The EHH test application defines a core region (eg, an allele with hypothesis of selection) and specifies an upstream and downstream distance from the variant. Then, calculates the probability of obtaining an identical core region by randomly taking two chromosomes [3].

Starting from the EHH statistic, Voight et al. developed the Integrated Haplotype Score (iHS) analysis in 2006, which compares the area under the curve defined by the EHH test [17]. Another variation of the EHH is the Cross Population Extended Haplotype Homozygosity (XP-EHH) statistic [18], which compares the length of haplotypes between populations. Other statistics based on linkage disequilibrium include the long-range haplotype (LRH)[19], LD decay (LDD)[20], and identity-by-descent (IBD) tests [21,22]. The focus of this work will be entirely directed at PBS and iHS.

Description of Population Branch Statistic (PBS)

Although the FST can be used to detect potential targets of natural selection, a simple ranking of FST values cannot reveal which of the two populations was affected by selection [11]. To overcome this issue, Yi and collaborators described for the first time the PBS.

This statistic includes two related populations - pop1 and pop2 - and a third, more distantly related population - popout -. Then, by comparing the three pairwise FST values between these three samples - pop1vspop2, pop1vspopout, and pop2vspopout - it can estimate the frequency change that occurred in the population of interest (pop1) since its divergence from its related population (pop2). Therefore, a population’s PBS value represents the amount of allele frequency change at a given locus since its divergence from the other two populations [11].

This approach is pretty similar to the LSBL metric, in which the FST is pair-wise calculated from three or more populations to identify population-specific changes in allele frequency given a genetic and demographic context [3]. Even though LSBL and PBS share the same basis, the last one uses the classical log transformation by Cavalli-Sforza [11,17].

Description of integrated Haplotype Score (iHS)

The integrated Haplotype Score (iHS) statistic is based on the Extended Haplotype Homozygosity (EHH) test first proposed in 2002 by Sabeti et al [5]. The EHH test estimates the probability that two randomly selected chromosomes are homozygous for all the variants found in a genetic distance X. The value of EHH is computed on a scale from 0 (no homozygosity, all haplotypes are different) to 1 (complete homozygosity, all haplotypes are the same). When an allele under selection increases its frequency, high levels of haplotypic homozygosity are observed. Therefore, in EHH versus Distance plots the area under the EHH curve will usually be larger for an allele under selection compared to that observed for a neutral allele.

Voight et al captured the previously described effect by calculating the integral of the observed decrease in homozygosity from a specific core region until reaching an EHH value equal to 0.05. This integrated EHH (iHH) value is denoted as iHHA if it is computed with respect to the ancestral allele and as iHHD if it is computed with respect to the derived allele. Thus, the iHS is obtained as follow: In(iHHA/iHHD).

Bioinformatics workflow: nf-selection

To automate the process of selection signatures detection with PBS and iHS, we developed nf-selection, a bioinformatics Nextflow pipeline that achieves the implementation of the analyses required to obtain results. The pipeline, and scripts developed for operation are publicly available via GitHub.

mkdir 1000G_data
cd 1000G_data
wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000_genomes_project/release/20190312_biallelic_SNV_and_INDEL/ALL.chr2.shapeit2_integrated_snvindels_v2a_27022019.GRCh38.phased.vcf.gz

wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call_samples_v3.20130502.ALL.panel

grep "CEU" integrated_call_samples_v3.20130502.ALL.panel | cut -f1 > ceu_individuals.txt
grep "CHB" integrated_call_samples_v3.20130502.ALL.panel | cut -f1 > chb_individuals.txt
grep "YRI" integrated_call_samples_v3.20130502.ALL.panel | cut -f1 > yri_individuals.txt

## We will create an additional file with all the samples
cat ceu_individuals.txt chb_individuals.txt yri_individuals.txt > all_samples.txt

vcftools --gzvcf out.vcf.gz --out filtered_vcf --keep all_samples.txt --chr 2 --from-bp 100000000 --to-bp 150000000 --recode

vcftools --vcf filtered_vcf.recode.vcf --out ceu_filtered --keep ceu_individuals.txt --recode

git clone https://github.com/fernanda-miron/nf-selection.git 
mv 1000G_data/filtered_vcf.recode.vcf nf-selection/test/data/pbs_files
mv 1000G_data/*.txt nf-selection/test/data/pbs_files
mv 1000G_data/ceu_filtered.recode.vcf nf-selection/test/data/ihs_files

cd nf-selection/test/data/ihs_files
nano design_file.csv

chromosome,path_vcf,path_genetic_map,path_ancestral,path_manifest
2,test/data/ihs_files/ceu_filtered.recode.vcf,test/data/ihs_files/chr2.b38.predicted.map,nf_modules/ancestral_fasta/ANCESTOR_for_chromosome_GRCh38_2_1_242193529_1.fa,nf_modules/manifest_annotation/manifest2

cd nf-selection/test/data/pbs_files
nano design_file.csv

path_vcf,path_pop1,path_pop2,path_popout
test/data/pbs_files/filtered_vcf.recode.vcf,test/data/pbs_files/ceu_individuals.txt,test/data/pbs_files/chb_individuals.txt,test/data/pbs_files/yri_individuals.txt

bash runtest.sh